///*
// * BatchReactor.cpp
// *
// *  Created on: 3 Jun 2011
// *      Author: allan
// */
//
//#include "BatchReactor.h"
//
//// C++ includes
//#include <cmath>
//#include <functional>
//#include <iostream> // TODO remove
//#include <iomanip> // TODO remove
//using namespace std::placeholders;
//
//// GSL includes
//#include <gsl/gsl_odeiv.h>
//
//// FluidReact includes
//#include "MultirootNewtonRaphson.h"
//#include "Speciation.h"
//#include "SystemOutput.h"
//
//BatchReactor::BatchReactor() :
//system(0)
//{}
//
//BatchReactor::~BatchReactor()
//{}
//
//void BatchReactor::SetSystem(System* system)
//{
//	this->system = system;
//}
//
//const System*const BatchReactor::GetSystem() const
//{
//	return system;
//}
//
//const vector<double> TotalConcentrationOfPrimarySpecies(const System* system)
//{
//	const unsigned numOfPrimarySpecies = system->GetNumOfPrimarySpecies();
//	
//	const unsigned numOfEquilibriumSpecies = system->GetNumOfEquilibriumSpecies();
//	
//	vector<double> Uj(numOfPrimarySpecies);
//	
//	for(unsigned iPrimary = 0; iPrimary < numOfPrimarySpecies; ++iPrimary)
//		Uj[iPrimary] = system->GetMolesOfSpecies(system->GetIndexOfPrimarySpecies(iPrimary));
//	
//	for(unsigned iEquilibrium = 0; iEquilibrium < numOfEquilibriumSpecies; ++iEquilibrium)
//	{
//		for(unsigned iProduct = 0; iProduct < system->GetEquilibriumReaction(iEquilibrium).GetNumOfProducts(); ++iProduct)
//		{
//			const unsigned j = system->GetEquilibriumReaction(iEquilibrium).GetIndexOfProduct(iProduct);
//			
//			const double vej = system->GetEquilibriumReaction(iEquilibrium).GetStoichiometryOfProduct(iProduct);
//			
//			Uj[j] += vej * system->GetMolesOfSpecies(system->GetIndexOfEquilibriumSpecies(iEquilibrium));
//		}
//	}
//	
//	return Uj;
//}
//
//struct BatchReactionParameters
//{
//	const System* system;
//
//	const vector<double>& ratesOfKineticSpecies;
//
//	const vector<double>& ratesOfPrimarySpecies;
//};
//
//int ODESystemFunction(double t, const double y[], double dydt[], void* params)
//{
//	BatchReactionParameters p = *(BatchReactionParameters*)params;
//
//	const unsigned numOfPrimarySpecies = p.system->GetNumOfPrimarySpecies();
//
//	const unsigned numOfKineticSpecies = p.system->GetNumOfKineticSpecies();
//
//	for(unsigned j = 0; j < numOfPrimarySpecies; ++j)
//	{
//		dydt[j] = p.ratesOfPrimarySpecies[j];
//	}
//	
//	for(unsigned k = 0; k < numOfKineticSpecies; ++k)
//	{
//		dydt[k + numOfPrimarySpecies] = p.ratesOfKineticSpecies[k];
//	}
//	
//	return GSL_SUCCESS;
//}
//
//void BatchReactor::Solve(double tEnd, double h, double hMax, string filename)
//{
//	const unsigned numOfPrimarySpecies = system->GetNumOfPrimarySpecies();
//	
//	const unsigned numOfKineticSpecies = system->GetNumOfKineticSpecies();
//
//	//const gsl_odeiv_step_type * type = gsl_odeiv_step_rk8pd;
//	const gsl_odeiv_step_type * type = gsl_odeiv_step_rk4;
//
//	gsl_odeiv_step* s    = gsl_odeiv_step_alloc(type, numOfPrimarySpecies + numOfKineticSpecies);
//	//gsl_odeiv_control* c = gsl_odeiv_control_y_new(0.0, 1.0E-10);
//	gsl_odeiv_control* c = 0;gsl_odeiv_control_yp_new(1.0E-20, 1.0E-20); // limiting derivatives (todo you need to understand the adaptive algorithm better)
//	gsl_odeiv_evolve* e  = gsl_odeiv_evolve_alloc(numOfPrimarySpecies + numOfKineticSpecies);
//
//	gsl_odeiv_system sys = {ODESystemFunction, 0, numOfPrimarySpecies + numOfKineticSpecies, 0};
//
//	double t = 0.0;
//	double* y    = new double[numOfPrimarySpecies + numOfKineticSpecies];
//	double* dydt = new double[numOfPrimarySpecies + numOfKineticSpecies];
//
//	vector<double> Ck(numOfKineticSpecies);
//	
//	for(unsigned iKinetic = 0; iKinetic < numOfKineticSpecies; ++iKinetic)
//	{
//		Ck[iKinetic] = system->GetMolesOfSpecies(system->GetIndexOfKineticSpecies(iKinetic));
//	}
//	
//	vector<double> Uj = TotalConcentrationOfPrimarySpecies(system);
//
//	for(unsigned j = 0; j < numOfPrimarySpecies; ++j)
//	{
//		y[j] = Uj[j];
//	}
//
//	for(unsigned k = 0; k < numOfKineticSpecies; ++k)
//	{
//		y[k + numOfPrimarySpecies] = Ck[k];
//	}
//	
//	SystemOutput sysOutput(system, "s");
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(*system);
//	
//	for(unsigned iPrimary = 0; iPrimary < numOfPrimarySpecies; ++iPrimary)
//		speciation.ImposeUj(iPrimary, Uj[iPrimary]);
//	
//	while(t < tEnd)
//	{
//		cout << fixed;
//		cout << "Progress: " << t/tEnd * 100 << "%" << "\t";
//		
//		sysOutput.Register(t, system);
//		
//		vector<double> ratesOfKineticSpecies = system->GenerationRatesOfKineticSpecies();
//		
//		vector<double> ratesOfPrimarySpecies = system->GenerationRatesOfPrimarySpecies(ratesOfKineticSpecies);
//		
//		BatchReactionParameters parameters = {system, ratesOfKineticSpecies, ratesOfPrimarySpecies};
//
//		sys.params = &parameters;
//		
//		int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &t, tEnd, &h, y);
//		
//		if(h > hMax) h = hMax;
//		
//		for(unsigned j = 0; j < numOfPrimarySpecies; ++j)
//		{
//			Uj[j] = y[j];
//		}
//
//		for(unsigned k = 0; k < numOfKineticSpecies; ++k)
//		{
//			Ck[k] = y[k + numOfPrimarySpecies];
//		}
//		
//		system->SetMolesOfKineticSpecies(Ck);
//		
//		for(unsigned iPrimary = 0; iPrimary < numOfPrimarySpecies; ++iPrimary)
//			speciation.ImposeHj(iPrimary, Uj[iPrimary]);
//		
//		bool converged = speciation.Solve(1.0E-6, 100);
//		
//		if(!converged)
//		{
//			sysOutput.Print(filename);
//			cout << "Error: The numerical solution of the nonlinear system did not converge." << endl;
//			exit(0);
//		}
//		
//		if(status != GSL_SUCCESS)
//		{
//			cout << "Failure in GSL ODE solver. Check what happened, Allan." << endl;
//			break;
//		}
//		
//		cout << endl;
//	}
//	
//	gsl_odeiv_evolve_free(e);
//	gsl_odeiv_control_free(c);
//	gsl_odeiv_step_free(s);
//
//	delete [] y;
//	delete [] dydt;
//	
//	sysOutput.Print(filename);
//}
